### Establish the out directory

out.dir = '/home1/99/jc152199/ChapterOne/DailyTmaxDistributionWithPhysCTMax/'

### List the files with climate extract from individual days

clim.dir = '/home1/99/jc152199/COPclimateextracts/'
clim.files = list.files(clim.dir,pattern='.csv',full.names=T)

### Get a list of unique species, removing generic 'COPSPP' records

a = read.csv(clim.files[1], header=T)

spp.list = unique(a$spp)[-2]

### Blank object to bind data onto

outdata = NULL

### Now begin looping through each clim.file, extracting max temp for each species from that day

for (ii in c(1:length(clim.files)))

	{
	
	### Read in a climate file
	
	tdata = read.csv(clim.files[ii], header=T)
	
	### Remove records pertaining to generic Cophixalus
	
	tdata = tdata[which(tdata$spp!='COPSPP'),]
	
	### Identify the hottest temperature experienced by each species according to AWAP and BRT
	
	BRT.maxdata = aggregate(tdata$BRT_tmax, by=list(spp=tdata$spp), FUN=max)
	AWAP.maxdata = aggregate(tdata$AWAP_tmax, by=list(spp=tdata$spp), FUN=max)
	
	### Merge two dataframes together and change some names
	
	maxdata = merge(BRT.maxdata, AWAP.maxdata, by=c('spp'))
	names(maxdata)[2] = 'BRT_tmax'
	names(maxdata)[3] = 'AWAP_tmax'
	
	### Bind data for writing out
	
	outdata = rbind(outdata,maxdata)
	
	cat('\n',clim.files[ii],'\n')
	
	}
	
### Write out outdata

write.csv(outdata, file=paste(clim.dir,'DailyTmaxForAllCophixalus.csv',sep=''), row.names=F)

#### Read in outdata

outdata = read.csv(paste(clim.dir,'DailyTmaxForAllCophixalus.csv',sep=''),header=T)

### Remove NA's

outdata = na.omit(outdata)

### Read in raw physiology data from base.dir

minmax = read.csv('/home1/99/jc152199/MAXENT/output/PhysData/minmax.csv', header=T)

#### Now need to loop through each species, plotting a density plot of outdata
#### Then overlay lines for CI's
#### Then overlay lines for CTmin and CTmax


for (spp in c('COPAENI','COPBOMB','COPCONC','COPNEGL','COPINFA','COPHOSM','COPEXIG','COPMONT'))

	{
	
	### First obtain the phys results for the species of interest
	
	spp.mm = minmax[which(minmax$Species==spp),]
	
	### Now subset outdata to only the species of interest
	
	spp.data = outdata[which(outdata$spp==spp),]
	
	### Calculate density objects
	
	BRT.density = density(spp.data[,2])
	AWAP.density = density(spp.data[,3])
	
	### Calculate the y lims
	
	ylims = c(0,max(c(BRT.density$y,AWAP.density$y)))
	xlims = range(c(BRT.density$x,AWAP.density$x))
	
	### Calculate a polygon that represents CTMax 95% CI (based on normally distributed data)
			
			if(nrow(na.omit(spp.mm))!=0)
			
			{
			
			max.poly = data.frame(x=c(spp.mm[1,5]-spp.mm[1,6]*1.96,spp.mm[1,5]-spp.mm[1,6]*1.96,spp.mm[1,5]+spp.mm[1,6]*1.96,spp.mm[1,5]+spp.mm[1,6]*1.96),y=c(ylims[2],0,0,ylims[2]))
			
			### Establish limits based on bc data and mc data
	
			xlims = range(c(BRT.density$x,AWAP.density$x,max.poly[,1]))
			
			}
			
	
	#### Now open the .png device driver and create a file for each species (but all in the same folder)

	png(paste(out.dir,spp,'2.png',sep=''),units='cm', height=18, width=18, res=600)
	
	### Configure the plot space, and plot the BRT density object
	
	plot(BRT.density, col='blue', xlim=c(xlims[1]-6,xlims[2]), ylim=ylims, sub='', main=paste(spp),xlab='Daily Max Temp', ylab='Density')
	
	### Add the AWAP density object
	
	points(AWAP.density, type='l', col='red')
	
	### Add line CTMax
	
	points(c(spp.mm[1,5],spp.mm[1,5]),c(0,ylims[2]),type='l',lty=3, col='#FF0000')
	
	### Add a legend
	
	legend('topleft',legend = c('Empirical Tmax Distribution','AWAP Tmax Distribution','CTMax 95% CI'),text.col=c('blue','red','black'),lty=c(1,1,1),col=c('blue','red','#FF000025'),lwd=c(1,1,8), bty='n')
	
	### Plot the polygon if it exists
	
			if(nrow(na.omit(spp.mm))!=0)
			
			{
			
			polygon(x=max.poly[,1],y=max.poly[,2], col='#FF000025', border=NA)
			
			}
	
	### Close the device driver
	
	dev.off()
	
	}
	
# Close loop
	
	
	
